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ABSTRACT 


Mean  gravity  anomalies  are  predicted  to  provide  the  average  value  of 
gravity  within  1°  x  1°  surface  areas  or  surface  areas  of  other 
sizes.  Prediction  methods  include  conventional,  statistical,  deter¬ 
ministic,  and  geophysical.  The  statistical  methods  appear  to  give  the 
best  results  when  adequate  computer  resources  are  available.  The 
deterministic  methods  give  good  results  with  minimum  time  and  computer 
requirements.  The  geophysical  methods  give  usable  predictions  when 
little  or  no  measured  data  is  available.  The  most  accurate  conven¬ 
tional  method  is  graphical  interpolation,  but  this  method  is  too  time- 
consuming  to  be  recommended  for  general  use. 

DEFINITION  AND  SCOPE 


Gravity  prediction  is  defined  as  any  process  that  enables  the  estima¬ 
tion  of  a  gravity  anomaly  value  (1)  for  a  site  where  gravity  has  not 
been  measured  or  (2)  that  represents  the  average  value  of  gravity 
within  a  given  surface  area. 

At  the  Defense  Mapping  Agency  Aerospace  Center  (DMAAC),  gravity 
prediction  is  used  almost  exclusively  for  estimation  of  mean  gravity 
anomalies,  that  is,  for  estimation  of  gravity  anomaly  values  that 
represent  the  average  values  of  gravity  within  a  given  surface  area. 
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Mean  gravity  anomalies  are  predicted  for  1*  x  1 ' ,  5'  x  5',  15'  x  15', 
30*  x  30'  ,  1°  x  1°/  and  5°  x  5°  surface  areas.  In  order  to  limit  the 
scope  of  this  paper,  we  will  confine  the  discussion  to  1°  x  1°  mean 
gravity  anomaly  prediction.  However,  most  of  the  methods  discussed 
are  also  suitable  for  prediction  of  values  that  represent  surface 
areas  of  other  sizes,  and  many  can  be  used  for  prediction  of  gravity 
at  specific  sites. 

The  discussions  of  gravity  prediction  in  this  paper  are  limited  to 
land  areas.  Different  results  may  be  obtained  when  these  methods  are 
applied  in  oceanic  areas. 

VALUE  OF  1°  X  1°  MEAN  GRAVITY  ANOMALIES 

A  continuous  gravity  anomaly  field  covering  the  earth* s  surface  is 
essential  for  solving  many  of  the  problems  of  physical  geodesy.  For 
example,  with  a  continuous  gravity  anomaly  field,  geoid  undulations 
and  deflections  of  the  vertical  can  be  determined  at  the  earth's 
surface,  and  gravity  disturbance  components  can  be  determined  in  space 
above  the  earth’s  surface. 

The  classical  integral  formulas  of  physical  geodesy,  such  as  those  of 
Stokes  and  Vening  Meinesz,  must  be  evaluated  by  summation.  In  these 
formulas,  small  finite  compartments  replace  infinitesimal  surface 
elements  in  the  summation.  These  compartments  are  obtained  by  sub¬ 
dividing  the  surface  of  the  earth  by  concentric  circles  and  their 
radii,  or  by  grid  lines  along  the  parallels  and  meridians  to  form 
"square"  blocks,  for  example,  1°  x  1*  blocks.  The  "square"  blocks  are 
summed  for  a  major  portion  of  the  earth's  surface.  The  1°  x  1°  blocks 
are  a  most  important  element  in  these  summations. 

In  addition,  a  worldwide  set  of  1°  x  1°  mean  gravity  anomalies  can  be 
used  to  represent  the  potential  of  the  earth.  Alternatively,  the  1°  x 
1*  mean  gravity  anomalies  can  be  summed  in  a  suitable  algorithm  to 
obtain  a  set  of  geopotential  harmonic  coefficients  that  represent  the 
earth’s  potential. 

The  1p  x  1°  mean  gravity  anomaly  values  needed  for  the  above  applica¬ 
tions  all  must  be  obtained  by  some  suitable  gravity  prediction 
ods. 

GRAVITY  PREDICTION  METHODS- INTRODUCTION 

In  general,  a  1°  x  1°  mean  gravity  anomaly  is  predicted  as  a  function 
of  gravity  that  is  measured  at  discrete  points  on  the  surface  of  the 
earth.  If  all  of  the  measured  gravity  data  lies  within  the  1°  x  1° 
area  for  which  the  prediction  is  required,  the  prediction  process  is 
essentially  interpolation  or,  in  some  cases,  merely  data  averaging. 
If  most  of  the  measured  data  lies  outside  of  the  1°  x  1°  area  in 
question,  the  prediction  becomes  essentially  an  extrapolation 
process.  Both  interpolation  and  extrapolation,  as  well  as  data 
averaging,  are  included  in  many  gravity  prediction  methods. 

There  is  really  only  one  way  to  obtain  a  rigorously  correct  1°  xl® 
mean  gravity  anomaly.  To  do  this,  it  is  necessary  to  apply  a  data 
averaging  integral  of  the  form  {Heiskanen  and  Moritz,  1967): 


rr  /  f  Ag(x,y)  dx  dy 
x-0  y=0 


ab 


(1) 


where 


Ag  =  1°  x  1°  mean  gravity  anomaly 

a,b  =  rectangular  dimensions  of  the  1°  x  1°  surface  are 

for  which  Ag  is  required, 

Ag  =  gravity  anomaly  values  computed  from  measured 

gravity  at  discrete  points  (x,y)  within  the  1°  x  1°  area. 

If  the  Ag  are  Bouguer  gravity  anomalies,  the  Ag  is  a  1°  x  1  °  mean 
Bouguer  anomaly.  If  the  Ag  are  free-air  gravity  anomalies,  the  Ag  is 
a  1°  x  1°  mean  free-air  gravity  anomaly.  Either  anomaly  form  can  be 
used  successfully  in  (1). 

In  order  to  apply  the  integral  (1),  the  measured  gravity  values  must 

be  available  for  all  points  (x,y),  differentially  small  distances 

apart,  throughout  the  1°  x  1°  area.  Of  course,  this  condition  is 
never  fulfilled  in  practice,  and  the  integral  (1)  can  never  be  applied 
as  written. 

The  automated  gravity  data  holdings  of  the  DoD  Gravity  Library  (DoDGL) 
exceed  9.5  million  measurements.  If  these  measurements  were  to  be 
divided  equally  among  the  64,800  1°  x  1°  areas  of  the  world,  there 

would  be  slightly  more  than  145  measurements  per  1°  x  1°  area. 
However,  the  9.5  million  measurements  actually  are  very  unequally 
divided  among  the  1°  x  1°  areas.  The  number  of  stations  per  1°  x  1p 
area  varies  from  a  maximum  of  19,400  to  a  minimum  of  zero.  Neverthe¬ 
less,  a  1°  x  1°  land  area  containing  more  than  100  measurements  must 
be  considered  unusual  (except  within  the  U.S.  and  parts  of  Western 
Europe).  In  fact,  most  land  areas  contain  a  few  tens  of  stations  or 
less.  Moreover,  the  measurements  within  a  1°  x  1°  area  often  are 
poorly  distributed.  Thus,  the  problem  of  gravity  prediction  usually 
is  to  obtain  the  best  estimate  of  the  mean  anomaly  based  upon  a 
limited  amount  of  unequally  distributed  measured  data. 

GRAVITY  PREDICTION  METHODS -DESCRIPTION 

Gravity  prediction  methods  that  have  been  used  generally  fall  into 
four  categories:  conventional,  statistical,  geophysical,  and  deter¬ 

ministic.  The  conventional  and  geophysical  methods  are  useful  only 
for  mean  anomaly  predictions.  The  statistical  and  deterministic 
methods  can  be  used  either  for  mean  anomaly  prediction  or  for  predic¬ 
tion  of  gravity  at  sites  where  gravity  has  not  been  measured. 

Graphical  Interpolation 

The  graphical  interpolation  method,  one  of  the  oldest  of  the  conven¬ 
tional  methods,  requires  that  a  gravity  anomaly  map  be  constructed 
from  the  measured  data.  The  gravity  anomaly  map  may  be  constructed 
either  by  automated  plotter  or  by  hand. 

A  gravity  anomaly  map  drawn  by  automated  plotter  is  considerably  more 
economical  than  such  a  map  drawn  by  hand.  The  machine  contouring 
algorithm  can  be  based  either  on  linear  interpolation  between  data 
points  or  on  least  squares  prediction.  In  areas  containing  ample 
amounts  of  well  distributed  measured  data,  machine  contouring  cer¬ 
tainly  is  just  as  reliable  as  hand  contouring.  However,  machine 
contours  may  be  drawn  inaccurately  over  areas  containing  relatively 
small  amounts  of  measured  data  and/or  poorly  distributed  data. 


Gravity  anomaly  contour  maps  drawn  by  hand  often  provide  a  good 
representation  -for  data  averaging  purposes-  of  the  gravity  field  even 
when  measured  data  is  relatively  sparse  and/or  poorly  distributed.  In 
addition,  such  maps  frequently  can  be  improved  by  using  data  from 
topographic  and  geologic  maps.  When  such  supplementary  information  is 
used,  the  placement  of  gravity  contours  is  influenced  by  considering 
topographic  and  geologic  structure  effects  on  gravity  anomaly  varia¬ 
tions.  Similar  techniques  may  be  used  to  enhance  gravity  anomaly 
contour  maps  drawn  by  automatic  plotters.  Hand  contouring,  although 
effective,  is  a  very  time-consuming  and  manpower  intensive  process, 
and  the  reliability  of  the  resulting  map  depends  upon  the  experience 
of  the  compiler.  The  contour  interval  selected  and  the  amount  of 
smoothing  used  in  drawing  the  contours  also  influence  the  quality  of 
the  results. 

Once  the  gravity  anomaly  map  has  been  drawn,  the  integration  (1)  can 
be  performed  graphically  with  reference  to  the  map.  One  method  of 
graphical  integration  is  to  average  contours  by  eye  within,  say,  5*  x 
5'  blocks.  Then  the  144  5'  x  5*  averages  within  the  1°  x  1°  block  are 
themselves  averaged  to  obtain  the  final  1°  x  1°  mean  gravity  anomaly 
prediction.  An  alternate  method  is  to  read  center  values  for  each  5' 
x  5'  area  from  the  contour  map,  and  then  average  the  144  center 
readings.  Tests  (undocumented)  at  DMAAC  have  shown  essentially  no 
difference  in  results  between  the  two  averaging  methods.  Conse¬ 
quently,  the  latter  method  is  preferable  because  of  its  simplicity  and 
because  it  requires  less  subjective  judgment. 

Predictions  in  continental  areas  always  are  made  using  Bouguer  gravity 
anomalies  because  this  anomaly  form  is  nearly  independent  of  corre¬ 
lation  with  short  wavelength  elevation  variations  and,  therefore,  well 
suited  for  the  purposes  of  interpolation  and  extrapolation.  The  mean 
Bouguer  gravity  anomaly  prediction  must  be  converted  to  a  mean  free- 
air  gravity  anomaly  which  is  the  appropriate  anomaly  form  for  geodetic 
applications.  The  mean  free-air  gravity  anomaly  is  computed  using 


Ago  «  Ag_  +  0. 1 1 19h 
F  B 


where 


A  gF  *  1°  x  1°  mean  free-air  gravity  anomaly 


Ag0  31  1°  k  1°  mean  Bouguer  gravity  anomaly 

B 


h  =*  mean  elevation  of  the  1°  x  1°  area  in  meters 
Least  Squares  Prediction 

There  are  a  limited  number  of  gravity  measurements  scattered  through  a 
typical  1°  x  1°  area,  rather  than  the  infinite  number  required  by  the 
integral  ( 1 ) .  A  reasonable  approach  to  this  situation  is  to  approxi¬ 
mate  the  true  mean  gravity  anomaly  by  some  linear  combination  of  the 
measured  values  using 


Ag  -  \  a.  Ag, 

i-1  1 


(3) 


where  __ 

Ag  =  1°  x  1°  mean  gravity  anomaly 

Ag.  =  gravity  anomaly  value  computed  from  measured  data 
for  the  ith  point 

=  arbitrary  coefficient  to  be  determined. 

The  coefficients,  a which  depend  upon  the  relative  positions  of  the 
gravity  measurements  and  mean  anomaly  value,  can  be  obtained  in  a 
number  of  different  ways.  Perhaps  one  of  the  best  choices  is  to 
determine  the  coefficients  in  such  a  way  that  the  standard  error  of 
prediction  is  minimized.  Such  a  determination  leads  to  the  least 
squares  prediction  (or  least  squares  estimation)  method. 

Least  squares  prediction  has  been  the  most  widely  used  statistical 
prediction  method.  The  formulation  for  least  squares  prediction  was 
developed  by  Moritz  and  later  modified  by  Rapp  for  practical  applica¬ 
tion  on  digital  computers.  Details  can  be  found  in  Heiskanen  and 
Moritz  (1967),  and  Rapp  (1964). 


Although  least  squares  prediction  can  be  used  to  predict  mean  anoma¬ 
lies  directly,  most  have  used  a  form  that  predicts  the  value  of 
gravity  at  a  discrete  point.  In  the  latter  case,  Moritz  shows  that, 
for  least  standard  error  of  prediction,  the  coefficients  of  (3) 
assume  the  form 
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where 


Ag  =  predicted  gravity  anomaly  at  a  point,  P 
P 

(-1) 


'ik 


=  inverse  of  a  matrix  that  expresses  the  autocovar¬ 


iance  between  all  gravity  anomalies  in  the  vicinity  of  the  point,  P. 


C  =  row  matrix  that  expresses  the  autocovariance 
between  theF  gravity  anomaly  at  the  point  P  and  the  gravity  anomalies 
used  in  the  prediction. 


Agk  *  column  matrix  of  the  measured  gravity  anomalies 
used  in  the  prediction. 


Both  interpolation  and  extrapolation  may  be  involved  in  the  applica¬ 
tion  of  ( 5 ) . 
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the  point  P 


square  of  the  standard  error  of  prediction  at  the 
variance  of  the  gravity  anomalies  in  the  vicinity  of 


It  can  be  seen  from  equation  (5)  that  least  squares  prediction  does 
estimate  the  value  of  the  gravity  anomaly  at  a  point.  To  obtain  a  1° 
x  1°  mean  gravity  anomaly,  equation  (5)  is  used  to  predict  gravity 
anomaly  values  at  the  center  of  5*  x  5*  blocks.  Then  the  144  5'  x  5* 
predictions  within  the  1°  x  1°  block  are  averaged  to  obtain  the  final 
1°  x  1°  prediction.  It  should  be  noted  that  the  averaging  of  predic¬ 
tions  at  5'  x  5*  block  centers  is  identical  to  the  procedure  generally 
used  to  obtain  the  final  1°  x  1°  prediction  in  the  graphical  interpo¬ 
lation  method. 


Least  squares  prediction  depends  upon  the  availability  of  a  reasonably 
representative  gravity  anomaly  autocovariance  function.  The  autoco¬ 
variance  function  represents  the  average  rate  of  change  of  the  gravity 
anomaly  field  as  a  function  of  distance  within  a  limited  region,  say, 
5°  x  5°.  It  is  important  that  a  local  rather  than  a  global  covariance 
function  be  used  for  gravity  anomaly  prediction  purposes.  This  will 
insure  compatibility  between  the  gravity  anomaly  covariance  coeffi¬ 
cients  used  for  prediction  and  the  gravity  field  characteristics  of 
the  prediction  area. 

The  local  covariance  function  for  an  area  often  can  be  represented  by 
an  analytical  expression,  first  suggested  by  Hirvonen  (1962),  of  the 
form 

C 

C(s)  -  - 2—r  (7) 

1+(B/d) 

where 

C(s)  ■  local  auto'rovar lance  of  the  gravity  anomalies  for  the 
distance,  a 

Cq  ■  local  variance  of  the  gravity  anomalies 

d  *  arbitrary  fixed  distance  that  enables  equation  (7)  to 

fit  the  local  covariance  function. 


Alternatively,  the  local  covariance  function  can  be  expressed  by  a 
polynomial  fit  to  the  data,  as  Rapp  (1964)  mentions,  using  a  formula 
of  the  form 
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where 


C1  #  *  *  cn  -  coefficients  that  enable  (8)  to  represent  the 

local  covariance  function 


In  order  to  be  used  in  the  automated  program  (Rapp,  1964)  generally 
applied  for  least  squares  prediction,  the  point  gravity  anomalies  used 
for  the  prediction  must  be  separated  by  a  minimum  of  0.4  minutes  of 


This  is  done  to  prevent  singularity  in  matrix  inversion.  Known 
systematic  influences  are  eliminated  by  deletion  or  scaling. 

Least  squares  prediction  is  among  the  most  widely  used  of  the  gravity 
anomaly  prediction  methods  because  it  is  completely  automated  and 
yields  results  which  are  theoretically  "best."  However,  when  predict¬ 
ing  gravity  anomalies  in  areas  containing  large  data  sets,  a  large 
amount  of  computer  memory  is  required. 


Least  Squares  Collocation 

Least  squares  collocation  is  a  natural  extension  of  least  squares 
prediction.  Like  least  squares  prediction,  least  squares  collocation 
depends  upon  the  statistical  properties  of  the  gravity  anomaly 
signal.  In  its  purest  form,  however,  least  squares  collocation 
combines  linear  prediction  with  a  consideration  of  the  geodetic 
boundary  value  problem. 


Least  squares  collocation  is  of  particular  interest  in  prediction 
problems  because  it  enables  inclusion  of  various  types  of  data,  not 
just  gravity  anomalies,  in  the  prediction  algorithm.  In  theory,  at 
least,  a  prediction  can  be  based  on  not  only  measured  gravity  data  but 
also  topographic  data,  geologic  and  geophysical  data,  and  the  var¬ 
iances  (accuracies)  of  these  data. 


However,  because  of  the  difficulty  in  specifying  firm  functional 
relationships  between  gravity  anomalies  and  other  types  of  data,  or 
even  obtaining  enough  numerical  geological  and  geophysical  data  for 
use  in  a  prediction,  most  investigators  have  been  content  to  add  only 
the  variances  of  gravity  anomaly  data  to  the  least  squares  prediction 
algorithm.  Thus,  in  (5),  the  covariance  matrix  C ^  becomes  +  D^r 
where  is  a  diagonal  matrix  of  gravity  anomaly  error  variances. 
Then  the  gravity  anomaly  prediction  procedure  follows  exactly  that 
described  in  the  previous  section  for  least  squares  prediction. 

Least  squares  collocation  is  definitely  the  prediction  method  of  the 
future,  and  several  researchers  now  are  active  in  development  of  the 
method  for  practical  application. 


Simple  Average 

The  simplest  way  of  applying  the  basic  prediction  equation  (3)  is  to 


set  all  of  the  coefficients  equal  to  1/n.  This  gives  the  equation 
for  the  conventional  simple  average  method 


Ag 


T  Ag 

i-1 


(9) 


Equation  (9)  must  be  applied  with  respect  to  Bouguer  gravity  anoma¬ 
lies  •  If  free-air  anomalies  are  used,  a  small  correction  term  must  be 
added  to  account  for  the  dependence  of  free-air  anomalies  upon  local 
elevation  changes. 

Sparse  data  and/or  uneven  distribution  of  the  measured  data  causes  the 
mean  gravity  anomaly  predicted  by  the  simple  average  method  to  be 
strongly  biased  in  areas  of  high  freguency  gravity  variations.  With  a 
dense,  even  distribution  of  data,  strongly  negative  or  positive  areas 
within  the  1°  x  1°  blocks  should  not  unduly  affect  the  result. 

Modified  Average 

The  modified  average  method  has  been  used  at  DMAAC  for  rapid  computer 
estimation  of  1°  x  1°  mean  gravity  anomalies  on  a  world-wide  basis. 
In  this  method,  a  simple  average  mean  gravity  anomaly  is  computed 
using  (9)  for  each  10'  x  10*  area  that  contains  measured  gravity  data 
within  a  1°  x  1°  block.  These  10*  x  10*  mean  anomalies  are  then 
averc  jed  in  two  steps  to  obtain  the  final  1°  x  1°  mean  gravity  anomaly 
prediction. 

The  appropriate  equations  are 
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Ag  =  predicted  1Q  x  1°  mean  gravity  anomaly 

Ag^0  =  30*  x  30*  mean  gravity  anomaly 

Ag10  =  10*  x  10'  mean  gravity  anomaly 

n  =  number  of  measured  gravity  anomalies  in  the  10'  x 

10'  block 

I*  =  number  of  non-empty  10'  x  10'  blocks  in  the  30'  x 

30'  block 

J  =  number  of  non-empty  30'  x  30'  blocks  in  the  1°  x  1° 

block 

Again,  the  prediction  must  be  done  in  terms  of  Bouguer  gravity  anoma¬ 
lies  for  the  equations  to  be  applicable. 

The  modified  average  method  is  a  significant  improvement  over  the 
simple  average.  For  example,  if  a  1°  x  1°  block  has  dense  data  in  one 
corner,  but  little  data  elsewhere,  the  simple  average  may  be  seriously 
biased.  The  modified  average  method,  on  the  other  hand,  gives  equal 
weight  to  the  contribution  of  each  10'  x  10'  block  within  the  1*  x  1° 
area,  thereby  reducing  biasing  due  to  poorly  distributed  data. 

Satisfactory  results  are  often  obtained  from  the  modified  average 
method  provided  that  there  are  not  too  many  empty  10 1  x  10*  blocks. 
Indeed  unsatisfactory  results  can  be  obtained  when  data  is  too 
sparse.  The  chief  value  of  the  method  is  its  ease  and  simplicity  of 
formulation  which  enables  a  rapid  first  approximation  of  1#  x  1°  mean 
values. 


Deterministic  Methods 

There  are  a  large  number  of  deterministic  methods  that  might  be  used 
for  gravity  prediction.  In  all  such  methods,  a  smooth  bivariate 
function  is  constructed  that  takes  on  the  prescribed  values  at  the 
data  points.  Most  of  these  methods  amount  to  fitting  a  smooth  polyno¬ 
mial  surface  to  the  gravity  anomaly  data.  The  prediction  is  given  by 
the  value  of  the  surface  at,  say,  the  center  of  a  5'  x  5'  block. 

The  two  deterministic  methods  discussed  here,  the  bicubic  spline  and 
multi-quadratic  analysis,  are  the  ones  that  appear  to  be  best  suited 
for  gravity  anomaly  prediction. 

Bicubic  Spline 

The  bicubic  spline  approximation  is  basically  a  finite  element  method 
used  to  evaluate  geodetic  integral  formulas.  Functions  to  be  inte¬ 
grated,  such  as  gravity  anomalies,  are  represented  by  a  piecewise 
cubic  two-dimensional  polynomial.  The  polynomial  is  fitted  together 
mathematically  to  form  a  continuous  transition  from  one  element  to  the 
next.  This  method  combines  the  analytic  handling  of  polynomials  with 
a  smooth  representation  by  the  resulting  function. 

For  the  purpose  of  digital  computation,  the  earth* s  surface  is  divided 
into  blocks,  say  5°  x  5°,  bounded  by  latitudes  and  longitudes. 
Gravity  anomaly  data  within  the  area  of  interest  are  used  to  compute 
the  polynomial  coefficients  which  determine  the  bicubic  spline* 

One  good  program  (Michlik,  1977)  uses  a  5*  integration  grid  and  a 
minimum  curvature  routine  for  selection  of  data  points  having  a 
significant  change.  The  integration  grid  data  is  merged  with  the 
cubic  spline  (in  both  x  and  y)  to  predict  gravity  anomalies  at  the 
center  of  each  5*  x  5*  block.  Then  the  144  5*  x  5*  blocks  within  each 
1°  x  1°  area  are  averaged  to  obtain  the  tinai  1W  x  iw  mean  gravity 
anomaly  prediction. 

The  spline  is  represented  by  (Davis  and  Kontis,  1970) 

3  3 

S(x,y)  =  l  l  M^x-x  )m  (y-y.  ,)“  (13) 

where  m=0  n=0  J 

*  unknown  for  the  two-dimensional  bicubic  surface. 

Yj  =  Xj  =  subintervals 

The  bicubic  spline  handles  data  distribution  irregularities  well  since 
the  interpolation  surface  passes  through  all  the  data  points.  The 
results  obtained  from  this  method  usually  are  very  good  unless  the 
data  is  very  sparse,  and  the  method  is  quite  economical  from  the 
standpoint  of  computer  resources  required.  The  major  drawback  is  that 
no  prediction  error  function  is  available. 


Multiquadric  Analysis 

Multiquadric  Analysis  is  a  least  squares  prediction  method  which 
applies  the  hypothesis  that  smooth  or  arbitrary  mathematical  surfaces 
can  be  approximated  to  any  desired  degree  of  exactness  by  the  summa¬ 
tion  of  quadric  forms.  However,  the  multiquadric  functions  are  based 
upon  geometric  and  physJ  :ral  consic  '.ration  rather  than  stochastic 
processes  (Hardy,  1978). 


There  are  two  forms  of  multiquadric  functions  that  can  be  used  for 
gravity  anomaly  prediction,  non-harmonic  and  harmonic.  The  non¬ 
harmonic  form  fits  a  hyperboloid  to  the  data  points.  The  mathematics 
of  the  harmonic  form  depend  upon  the  harmonic  function  reciprocal 
distance.  The  harmonic  form  can  represent  physical  properties  such  as 
mass,  and  may  be  interpreted  in  terms  of  the  point  mass  summation 
concept. 

The  general  equations  are: 

Non-harmonic  form 
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Harmonic  form 
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where 


xj'yj 

xi'yi 


universal  gravitational  constant 
radius  of  the  outer  sphere 
radius  of  the  inner  sphere 
coordinates  of  the  nodal  points 
coordinates  of  the  data  points 
coefficients  to  be  determined 
co- latitude 
longitude 


The  coefficients  are  computed  from  all  known  gravity  anomaly  values 
within  an  interpolation  area  so  that  one  set  of  coefficients  applies 
for  the  whole  area.  The  derived  coefficients  then  are  used  to  predict 
gravity  anomaly  values  at  5*  x  51  intervals  that  can  be  averaged  to 
obtain  the  final  1°  x  1°  mean  gravity  anomaly  prediction. 
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A  constant  radius  (r)  for  a  so-called  inner  sphere  on  which  the  ooint 
masses  lie  is  required  for  computation  of  the  coefficients  as  well  as 
the  predicted  point  values.  A  table  of  recommended  r  values  for  n 
point  mass  anomalies  is  given  by  Hardy  (1978). 

Multiquadric  analysis  has  been  recommended  as  the  most  economical  and 
accurate  prediction  method  by  Hein  and  Lenze  (1979)  and  Franke 
(1979).  However,  it  should  be  noted  that  these  investigators  were 
predicting  elevations  rather  than  gravity  anomalies. 

Geophysical  Prediction 

In  areas  containing  very  little  measured  gravity  data,  extrapolation 
based  upon  correlations  between  gravity  anomaly  variations  and  the 
corresponding  variations  in  geological,  geophysical,  and  topographical 
parameters  can  be  developed.  Thus, 

6  Ag  =  f(6h,  8S)  (20) 

where 

6  Ag  =  variations  in  the  1°  x  1°  gravity  anomaly. 

f(<5h,  AS)  =  some  function  of  topographic  and  structural 
variation,  respectively. 

If,  for  example,  the  changes  in  the  regional  portion  of  the  1°  x  1° 
mean  gravity  anomalies  are  constant  with  respect  to  the  corresponding 
changes  in  regional  topography,  which  is  the  case  within  many  regions 
of  homogenious  geological  structure,  then 


6  h 

where  0  is  a  constant. 

Integration  of  (21)  gives  the  function 


Ag  =  Rh  +  a 


(22) 


where  a  is  also  a  constant. 

Equation  (22)  is  a  form  of  the  so-called  basic  prediction  that  is  used 
in  one  variety  of  geophysical  prediction.  A  number  of  corrections  are 
added  to  (22)  to  obtain  the  complete  prediction  formula. 


Ag  =  (ph  +  a  )  +  f^S)  +  f2(S)  +  f3(h)  (23) 


In  (23),  f i { S )  corrects  for  long  period  gravity  variation  that  violate 
the  constancy  assumption  inherent  in  (21).  The  corrections  f^fs)  an(1 
f^(h)  add  the  gravitational  effects  of  short  period  structural  and 
topographic  changes,  respectively,  that  affect  each  1°  x  1°  area  in  a 
different  way,  and,  therefore,  cannot  be  modeled  by  the  basic  predic¬ 
tor  (22). 
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Details  of  this  and  other  types  of  geophysical  prediction  methods  may 
be  found  in  Wilcox  (1974). 

The  geophysical  prediction  method  provides  useable  estimates  of  1°  x 
1°  mean  gravity  anomalies  in  areas  where  othe^  types  of  prediction 
methods  fail  due  to  inadequate  amounts  of  measured  data  being  avail¬ 
able.  However,  the  principles  of  geophysical  prediction  potentially 
are  valuable  to  improve  predictions  made  by  other  methods  even  when 
adequate  data  is  available.  The  most  promising  method  for  incorpora¬ 
tion  of  geophysical,  geological,  and  topographic  parameters  in  gravity 
prediction  is  least  squares  collocation. 

COMPARATIVE  ANALYSIS  OF 
SELECTED  PREDICTION  METHODS 

A  5°  x  5°  area  was  selected  in  the  United  States  to  test  five  differ¬ 
ent  prediction  methods.  The  criteria  used  for  selecting  the  area  were 
(1)  a  good  sample  of  elevation  ranges  and  (2)  an  abundance  of  well 
distributed,  dense  measured  gravity  anomaly  data.  The  area  found  that 
met  these  criteria  is  bounded  by  35°N-40°N  and  115°W-120°W. 

The  five  methods  selected  for  testing  include  (1)  the  modified  average 
(MA),  (2)  graphical  interpolation  (GI),  (3)  cubic  spline  (SP),  (4) 
simple  average  (SA),  and  (5)  least  squares  prediction  using  the 
algorithms  developed  by  Rapp  (RP).  We  had  originally  intended  to  test 
multiquadric  analysis  fully,  but  due  to  problems  arising  in  refining 
part  of  the  basic  formulation  coupled  with  time  constraints  it  was 
necessary  to  omit  this  method  from  the  test. 

The  test  was  done  in  two  parts.  In  part  1,  some  33,000  measured 
gravity  anomalies  having  a  relatively  even  distribution  were  used  in 
each  prediction  method  to  predict  1°  x  1°  mean  Bouguer  gravity  anoma¬ 
lies.  Where  required  by  the  method  (e.g.  cubic  spline),  Bouguer  point 
anomalies  were  first  predicted  at  the  centers  of  a  5'  x  5*  integration 
grid  within  each  1°  x  1°  area,  and  these  were  averaged  to  obtain  the 
final  1°  x  1°  mean  Bouguer  anomaly  prediction. 

In  part  2,  the  data  was  thinned  at  random  to  simulate  a  sparse  (moder¬ 
ate  to  extreme),  poorly  distributed  data  field.  A  total  of  6,600 
measured  gravity  anomalies  were  used.  The  mean  Bouguer  gravity 
anomalies  were  predicted  as  in  part  1. 

Using  the  predicted  1°  x  1®  mean  Bouguer  gravity  anomalies  and  1°  x  1° 
mean  elevatior  values,  25  1°  x  1°  mean  free-air  gravity  anomalies  were 
obtained  using  equation  (2)  for  each  data  set:  full  (part  1)  and 
sparse  (part  2). 

Least  squares  prediction  for  the  full  data  set  (part  1)  was  chosen  as 
the  standard  for  statistical  comparison.  The  1°  x  1°  mean  free-air 
gravity  anomalies,  which  were  predicted  by  the  various  methods  in  both 
part  1  and  part  2,  were  compared  with  those  predicted  by  least  squares 
prediction  in  part  1. 

Mean  free-air  gravity  anomalies  predicted  by  each  method  for  all  1°  x 
1®  areas  within  the  selected  5°  x  5°  area  are  tabulated  in  Table  1  for 
part  1  (full  density)  and  Table  2  for  part  2  (sparse  density).  The 
mean  difference  (x) ,  standard  deviation  (o),  and  maximum  difference 
between  predictions  made  by  each  method  and  Rapp  full  density  predic¬ 
tions  (used  as  the  standard)  are  given  in  Table  3. 


From  Table  3,  it  can  be  seen  that  least  squares  prediction  (Rapp), 
graphical  interpolation,  and  the  cubic  spline  agree  satisfactorily  in 
both  the  full  and  sparse  areas.  There  is  no  appreciable  difference  in 
the  mean  free-air  gravity  anomalies  predicted  by  these  three  methods 
in  part  1  (full  distribution),  and  the  differences  in  part  2  are 
relatively  small.  Thus  any  of  these  three  methods  might  have  been 
selected  as  the  standard  for  comparison.  It  should  be  noted  that,  in 
both  parts  1  and  2,  graphical  interpolation  compares  most  favorably 
with  least  squares  prediction. 

The  modified  average  method  compares  closely  with  least  squares 
prediction  in  part  1  but  not  in  the  sparse  case  (part  2).  Worst 
results  in  part  2  were  for  two  1°  x  1°  areas  having  exceptionally  poor 
distribution  and  extremely  sparse  data.  The  simple  average  does  not 
give  good  results  in  either  case,  full  or  sparse. 

Since  x  (Table  3)  is  uniformly  small,  there  does  not  appear  to  be  a 
serious  systematic  error  component  in  prediction  for  any  method. 

All  prediction  computations  were  made  in  a  UNIVAC  1100  series  comput¬ 
ing  system.  A  time  study  of  the  prediction  methods  tested  were  made 
(Table  4)  based  upon  the  following  criteria:  (1)  preparation  time 
(data  tape  preparation  and  computer  deck  set  up),  (2)  total  computer 
time  (CPU  +  I/O)  necessary  for  obtaining  results,  (3)  other  (time 
necessary  for  obtaining  results  other  than  by  computer)  -  graphical 
interpolation  was  the  only  method  having  time  requirements  in  this 
category,  and  (4)  total  time  -  sum  of  (1),  (2)  and  (3). 

The  prediction  method  requiring  the  least  amount  of  total  time  was 
judged  to  be  the  most  economical.  The  bicubic  spline  proved  to  be  the 
most  economical  in  both  the  full  and  sparse  areas.  The  other  predic¬ 
tion  methods  ranked  in  the  following  order,  (1)  modified  average/sim¬ 
ple  average,  (2)  least  squares  prediction  (Rapp),  and  (3)  graphical 
interpolation . 

It  is  unfortunate  that  multiquadric  analysis  could  not  be  included  in 
this  test.  Other  investigators  (Hein  and  Lenze,  1979  and  Franke, 
1979)  report  that  multiquadric  analysis  is  both  very  accurate  and  very 
economical . 

Based  on  our  own  test  and  those  of  others,  we  recommend  the  following: 

(1)  Where  adequate  computer  resources  and  time  are  available,  the 
most  advantageous  and  accurate  prediction  method  is  most  probably 
least  squares  prediction  or  its  close  cousin,  least  squares  colloca¬ 
tion. 


(2)  Where  computer  capacity  or  time  is  a  constraint,  either  of 
the  two  deterministic  methods  -  bicubic  spline  or  multiquadric  analy¬ 
sis,  should  give  good  results. 


(3)  Graphical  interpolation  gives  good  results  -  and  is  certainly 
one  of  the  most  accurate  methods  even  when  data  is  sparse  or  poorly 
distributed  -  but  cannot  be  recommended  for  general  use  because  it  is 
the  most  time  consuming,  by  far,  of  all  methods  tested. 


(4)  The  modified  average  method  should  be  used  with  caution.  It 
is  an  acceptable  method  to  obtain  a  rapid  first  approximation,  but 
larger  prediction  errors  may  occur  in  some  instances. 

(5)  The  simple  average  method  should  never  be  used. 

Table  1 

Predicted  Mean  Free-Air  Anomalies 
Part  1  (Full  Density) 


40  °N 


MA 

4.2805* 

-12.2089 

4.4102 

13.2902 

12.2522 

GI 

4.4597 

-12.9616 

4.3776 

13.8333 

12.3807 

SP 

3.7975 

-12.5089 

5.4142 

13.4632 

12.6642 

SA 

5.4805 

-13.1089 

4.2102 

12.5902 

11.5528 

RP 

3.9200 

-13.0600 

3.9300 

13.5800 

12.8100 

MA 

33.9871 

-  .7109 

8. 1029 

10.7806 

-  5.3134 

GI 

32.9580 

-  1.1269 

8.2730 

10.8584 

-  6.2669 

SP 

32.8531 

-  .8549 

8. 1469 

10.8048 

-  5.8544 

SA 

27.6871 

-  2.2109 

9.6029 

12.7806 

-  4.9134 

RP 

33.0300 

-  .9800 

8.2500 

10.9500 

-  6.0900 

MA 

22.2060 

33.8316 

-  7.1673 

-  1.9487 

-  7.2170 

GI 

20.7261 

34.2795 

-  6.7499 

-  1.6327 

-  6.2753 

SP 

20.9750 

34.1456 

-  6.9753 

-  1.1927 

-  6.3070 

SA 

32.4060 

31.0316 

-  8.5673 

-  .4487 

-  7.4170 

RP 

20.7600 

34.3000 

-  6.5200 

-  1.6400 

-  6.4200 

MA 

-35.1653 

47.4328 

-12.7870 

-21 .2781 

13.0389 

GI 

-35.5764 

46.8307 

-13.4363 

-22.4559 

13.4485 

SP 

-35.7333 

46.6088 

-13.0860 

-22.0861 

13.4829 

SA 

-39.2653 

42.4328 

-17.3870 

-14.3781 

15. 1389 

RP 

-35.5100 

46.9200 

-13.1300 

-21.4100 

13.4400 

MA 

-42.7607 

10.4637 

-17.3472 

-22.5675 

2.0923 

GI 

-42.4919 

10. 1630 

-16.7819 

-23.4418 

2.1972 

SP 

-42.7970 

9.9907 

-16.9043 

-23.1914 

1.9853 

SA 

-46.1607 

21 .36^7 

-23.3472 

-19.8675 

-  .5077 

RP 

-42.4300 

10.1300 

-16.9100 

-23.2300 

2.1100 

35°N 

120 

°W 

1 15°W 

MA 

Modified 

Average 

GI 

m 

Graphical 

Interpolation 

SP 

s 

Spline 

SA 

33 

Simple  Average 

RP 

xz 

Rapp 

*MGAL 


Table  2 


Predicted  Mean  Free-Air  Anomalies 
Part  2  (Sparse  Density) 


40  °N 


MA 

3.8756* 

-10.8977 

9.7163 

17.8529 

13.9890 

GI 

4.0395 

-11.8158 

7.1011 

16.7485 

12.0591 

SP 

4.0425 

-11.7919 

7.3122 

20.6552 

14.6592 

SA 

4.9692 

-15.2477 

8.2114 

17.9613 

14.1267 

RP 

3.9500 

-11.9500 

7.3100 

16.6600 

14. 1500 

MA 

34.8429 

-  1.2623 

7.2511 

7.8838 

-  4.0134 

GI 

36.3538 

3.1648 

6.1723 

12.4487 

-  .1870 

SP 

34.2551 

.2720 

5.9579 

12.0946 

-  1.1784 

SA 

32.1871 

2.6715 

11.7380 

5.0632 

-  6.8134 

RP 

34.9600 

-  .0100 

6.0700 

11 . 1400 

1.3600 

MA 

21.6139 

33.0059 

-10.4673 

1.0513 

-  4.8118 

GI 

22.1643 

36.0593 

-  9.5799 

3.5916 

-  5.6321 

SP 

21.2360 

34.2326 

-  9.8313 

5.1893 

-  4.6010 

SA 

23.6145 

29.5273 

-  8.3673 

14.1013 

-  5.8434 

RP 

20.8400 

34.8000 

-  8.7900 

2.8200 

-  5.3200 

MA 

-36.1111 

47.6720 

-33.1160 

-17.6604 

16.0125 

GI 

-35.5056 

47.0307 

-19.6301 

-19.8031 

14.6663 

SP 

-33.9279 

46.0418 

-16.3920 

-16.8341 

14.3569 

SA 

-32.5754 

39.9168 

-45.0033 

-18.9781 

19.0337 

RP 

-34.9500 

47.1800 

-19. 1600 

-15.1300 

14.0000 

MA 

-41.4611 

14.6848 

-28.7605 

-17.7605 

2.8341 

GI 

-41.7183 

11.7366 

-22.4715 

-22.5474 

5. 1867 

SP 

-44.6513 

12.5447 

-25.8952 

-22.7094 

3.2273 

SA 

-36. 1213 

28.2094 

-36.1972 

-15.8875 

7.3038 

RP 

-41.2200 

11.2600 

-22.5000 

-21.9600 

4. 1200 

35°N 

120°W  1 1 5°W 

MA  =  Modified  Average 

GI  =  Graphical  Interpolation 

SP  «  Spline 

SA  =*  Simple  Average 

RP  *  Rapp 


*MGAL 


Table  3 


Statistical  Analysis 

I 


Prediction  Method 

Density 
of  Data 

X 

(mgal) 

0 

(mqal ) 

Maximum 

Difference 

Modified  Average 

Full 

Sparse 

*1788 

-.0340 

.5528 

5.4674 

1.45 

19.99 

Graphical 

Full 

Sparse 

-.0489 

.9133 

.2971 

2.914 

1.05 

6.5 

Spline 

Full 

Sparse 

.0016 

.8598 

.4062 

3.3652 

1.48 

8.985 

Simple  Average 

Full 

Sparse 

.2753 

.4320 

4.4937 

9.8200 

11.65 

31  .87 

Rapp* 

Full 

Sparse 

* 

.9152 

* 

2.967 

6.28 

♦Used  as  the  standard 


Table  4 

Prediction  Methods  Time  Comparison 
Total 

Density  Preparation  Computer  Time 
Method  of  Data  Time  (CPU  +  I/O)  Other  Total  Time 

_ (hr,  min)  (hr,  min, sec) _ 


MA/SA 

Full 

0 

hr 

25' 

0 

hr 

5 ' 26 . 0" 

0 

hr 

30 1 26.0” 

Sparse 

0 

hr 

15' 

0 

hr 

2 ' 10 . 0" 

0 

hr 

17' 10.0” 

GI 

Full 

6 

hr 

O' 

0 

h 

2'  5.0" 

64 

hr 

2'* 

70 

hr 

4*  5.0” 

Sparse 

3 

hr 

O' 

0 

hr 

1*  2.0" 

24 

hr 

59  *  * 

28 

hr 

1*  2.0” 

SP 

Full 

0 

hr 

35  * 

0 

hr 

2' 10.780" 

0 

hr 

37*10.78” 

Sparse 

0 

hr 

35' 

0 

hr 

1*56. 6" 

0 

hr 

36 '56.5" 

RP 

Full 

1 

hr 

ID 

8 

hr 

20' 11.987" 

9 

hr 

55*11.978 

Sparse 

1 

hr 

20' 

3 

hr 

13'38.002" 

4 

hr 

33*38.002 

♦Time  used  to  predict  5*  x  5*  center  values  per  1°  square 
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